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1.  Introduction 


Suppose  that  y(s), s  >  0,  is  a  continuously  differentiable  process.  Let  u(s)  and  u(s),  u(s)  >  o{s). 
be  continuously  differentiable  barriers.  Assume  that  the  process  starts  between  the  barriers,  i.o. 
u(0)  <  2/(0)  <  «(0)  and  if  2/(0)  -  u( 0)  or  y(O)  =  u(0)  then  y'(0)  <  «'(()),  y'(0)  >  V'(Q),  respectively. 

In  this  paper,  we  are  interested  in  the  densities  of  the  absorption  times  Tu,  Tv  of  the  y-process  in 
the  barriers  u,v,  respectively.  More  precisely,  let  71',  T'  be  the  first  passage  times  to  the  barriers 
u,  v ,  respectively.  Then  the  absorption  times  Tu,  Tv  are  defined  as  follows 


'  rpl 

J-u 

< 

if  T'  <  T', 

,  +oo 

otherwise, 

{  K 

if  rv<ru, 

\  +oo 

otherwise. 

Our  main  result,  presented  in  Section  2,  states  that  the  densities  of  TU,TV  can  be  expressed  in 
terms  of  conditional  expectations  in  the  following  way 

(1)  /T.(0  =E  [/(o,t)(2/)(2/'(0  ~  «'(0)+  1 2/(0  =  «(<)]  /y(t)MO). 

/t„(0  -E  [/(o,t)(2/)(2/,(0  -  At)Y  1 2/(0  =  KO]  fy(t)Kt)), 
where  z+  =  max(0,x),  x~  =  max( 0, -x)  and  I(o,t)(y)  is  the  indicator  function  defined  equal  to  1 
if  the  sample  path  does  not  cross  the  barriers  u,  v  prior  to  time  t  and  equal  to  0  otherwise.  The 
formula  (1)  is  an  extension  of  Durbin’s  formula  for  the  first  passage  density  [3,  10],  which  can  be 
obtained  by  replacing  the  lower  barrier  v  by  -oo.  Since  the  indicator  I(o,t)  is  a  function  of  the  whole 
sample  path  of  the  y-process,  the  expectations  in  (1)  are  difficult  to  evaluate  exactly.  However,  we 
shall  use  (1),  in  Section  3,  to  construct  upper-  and  lower-bounds  for  the  densities  of  TU,T„. 

Several  important  applications  are  related  to  the  two  absorbing  barrier  problem.  Two  are  dis¬ 
cussed  in  Section  4.  The  first  one,  which  arises  in  oceanography,  is  the  evaluation  of  the  joint 
distribution  of  wave-length  and  amplitude  of  random  waves;  i.e.  the  joint  distribution  of  the  differ¬ 
ence  in  time  between  the  uperossing  of  the  mean  sea  level  and  the  following  downcrossing  of  this 
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level,  and  the  highest  value  of  the  sea  in  this  interval,  (see  Section  4.2).  In  the  second  application. 

discussed  in  Section  4.3,  we  give  approximations  for  the  distribution  of  the  so  called  rainflow  cycle 

amplitude.  (The  Rainflow  method  was  developed  in  fatigue  analysis  to  describe  a  load  process.) 

2.  Basic  theorem 

We  begin  with  a  definition  of  a  class  of  processes  for  which  (1)  holds. 

Definition  1.  Let  y(s)  be  a  continuously  differentiable  process.  Assume  that  there  exists  a  k- 
variate  continuously  differentiable  process  A (s),  and  a  random  variable  X,  independent  of  A (s). 
with  bounded  and  continuous  density  function  such  that 

(2)  y(s)  =  F(s,X,A(s)),  s  >  0, 

where  F  is  a  continuously  differentiable  mapping  from  R+  x  Rx  Rk  into  R.  For  a  fixed  t  >  0,  the 
y-process  will  be  called  decomposable  at  t  if  there  exists  e  >  0  such  that  for  all  s,  |s  -  <|  <  e,  and 
all  z  £  Rk,  F(s,  -,z)  is  one-to-one.  Further,  if  the  y-process  is  decomposable  at  almost  all  t  then  y 
is  called  decomposable. 

The  class  of  decomposable  processes  is  quite  large  and  contains  for  example:  Gaussian  processes, 
functions  of  Gaussian  vector  processes,  Slepian  model  processes,  the  sum  of  a  Gaussian  process  and 
any  independent  continuously  differentiable  process,  etc.  An  example  of  a  class  of  processes  which 
are  not  decomposable  are  processes  which  are  “deterministic” on  some  interval,  e.g.  y(.s)  =  (/{$). 
g(s)  is  a  continuously  differentiable  function. 

We  turn  now  to  the  definition  of  the  particular  version  of  conditional  expectations  used  in  (1). 
Assume  that  y  is  a  decomposable  process  at  t.  Denote  by  ptlZ  an  inverse  mapping  to  F(t,-,z)  (2). 
i.e. 

(3)  F(t,pt,z{r),z)  =  r, 
where  r  £  R.  Let  yr  be  the  following  process 

(4)  Vr(s)  =  F(s,pliA{l){r),  A(s)). 
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Lemma  2.  Assume  that  the  y-process  is  decomposable  at  t.  If  h  is  a  non  negative  measurable 
functional  defined  on  y,  then  for  any  r  6  R 


(5)  E[h(y)\y(t)  =  r]/y(t)(r)  =  E[h(yr)  •  /(r|A(f))] , 
where  the  process  yr  is  defined  by  (4),  and  the  function  f(r |z)  is  given  by 

(6)  f(r\z)  =  )  dPtQ~r- 1  -fx  (Pt,z(r)) . 

The  function  p  is  defined  by  (3)  and  fx  is  the  density  of  X. 

Proof:  The  lemma  follows  from  Fubini’s  theorem. 

Obviously,  for  many  processes  y,  the  decomposition  (2)  is  not  unique.  In  that  case  one  can  choose 
the  decomposition  (2),  which  gives  the  most  convenient  expression  for  the  process  yr  (4).  For 
example,  if  y  is  a  zero-mean  Gaussian  process,  with  Var(y(t))  >  0,  the  most  natural  decomposition 
is 

(7)  y(s)  =  y(t)bt($)  +  At(s), 

where  At(s)  is  a  zero-mean  Gaussian  residual  process  independent  of  y(t)  with  covariance  function 
*'t(si,s2)  =  Cov(y(si),y(so)\y(t)),  and  bt(s)  =  Cov(y(s),y(t))/Var(y(t)).  Now,  by  (4),  the  process 
yr  is  defined  by  yr($)  =  r  •  bt(s)  +  A*(.s).  Since,  pt>z(r)  =  r,  the  formula  (5)  can  be  written  as 

E[h(y)\y(t)  =  r}fy[t){r )  =  E[h(yr)]fy{l)(r). 

In  general,  h(yr)  and  f(r\ A(t))  in  (5),  are  dependent  random  variables.  Wo  shall  illuminate  this 
by  using  a  different  decomposition  (7),  e.g. 

y(s)  =  y(0)bo(s)  +  A0(s). 

Suppose  b0(t)  ^  0,  then,  by  (3),  pt,:(r)  =  ff(f)  an^  fche  process  yr  is  given  by 

yr (^)  =  r  -f  A0(s). 
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Finally,  the  conditional  expectation  (5)  can  be  written  as 


E[h{y)\y(t)  =  r]/y(0(r)  =  E 


h{yr)m\fm 


(  r  —  Ap(i)  V 

V  m  )y 


Theorem  3.  Assume  that  the  process  y  is  decomposable  at  t.  If  i?[|j/'(i)||j/(t)  =  u(t)}  <  +cc. 
jB[|j/'(<)||y(t)  =  u(t)]  <  +oo,  see  (5),  then  the  densities  ofTu,TV)  are  finite  and  are  given  by  by  (1 ). 


Proof:  The  proof  is  similar  to  the  proof  of  Theorem  2  in  [10]. 

3.  Bounds  for  the  absorption  times  densities 


3.1  Introduction.  In  this  section  we  present  upper-  and  lower-bounds  for  the  density  of  Tu;  the 
density  of  Tv  can  be  treated  similarly. 

For  a  fixed  time  point  s,  denote  by  I(y,s)  the  indicator  function  defined  equal  to  1  if  u(s)  < 
y(s)  <  u(s )  and  equal  to  0  otherwise.  In  the  same  way,  for  a  vector  of  time  points  s  =  (si,. . . 

0  <  s,  <  t,  let  be  the  following  indicator 

n 

(8)  I{y;si,...,sn)  =  1[[l{y\Si). 

1=1 

Since,  for  any  nctor  s,  I{o,t){y)  <  I{v\  s),  an  upper  bound  for  the  density  of  Tu  can  be  obtained  by 
replacing  in  (1)  the  indicator  J( o,t)  by  I(y,  s),  i.e. 


(9)  fu(^n)  =  -  «'(<))+|y(0  =  «(0]/y(o(«(0)- 


However,  it  is  in  general  difficult  to  give  useful  lower  bounds  for  the  indicator  I(o,t){y)  in  ( 1 ),  and 
therefore  formula  (1)  is  not  useful  in  the  construction  of  lower  bounds  for  the  density  of  Tu.  Hence, 
as  in  [13],  we  prove,  in  Theorem  4,  a  second  formula  for  the  density  of  Tu,  which  can  be  used 
to  construct  lower  bounds.  The  approach  taken  in  [13]  leads  to  very  general  lower  bounds,  but 
the  numerical  effort  to  evaluate  these  lower  bounds  is  much  bigger  than  that  corresponding  to 
the  upper  bound.  By  some  further  restriction  on  the  residual  process  A  in  (1),  e.g.  when  A  is 
a  Gaussian  process,  we  can  construct  lower  bounds  of  the  same  complexity  as  the  upper  bounds 
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(9) .  In  addition,  by  (2),  the  process  y  is  a  function  of  the  process  A  and  y'(t )  is  a  function  of 
A(<),  A'(f),  and  hence,  in  order  to  evaluate  numerically  the  upper  bound  (9),  the  joint  density  of 
A(s),  A(t),  A'(f)  must  be  given  in  an  explicit  form.  Thus,  from  this  point  on,  we  assume  that  A 
is  a  zero-mean  continuously  differentiable  vector  valued  Gaussian  process. 

3.2  A  second  formula  for  the  density  of  Tu.  Let  A  be  a  zero-mean  continuously  differentiable 
vector  valued  Gaussian  process.  For  a  fixet  t  >  0,  consider  the  following  sequence  of  random  vectors 

(10)  <<°>=(A(0,A'(i)), 


£<n)  =(  A(t),  A'(.<),  A(5i A(5„)) , 

where  Si,  0  <  Si  <  t,  i  =  1, . . . ,  n,  are  random  times,  such  that  each  Si  is  a  function  of  alone, 

i.e.  if  =  zf*'"1),  then  Si(z(‘-1))  =  s,,  where  0  <  s;  <  t  is  a  fixed  time  point. 

In  the  following  we  use  the  decomposition  of  the  process  A  into  the  conditional  expectation  on 
and  the  residual  process,  i.e. 

(11)  A(s)  =  i;fAiS)iC{n)]  +  An(«). 

In  (11)  Lemma  9,  we  gave  an  explicit  formula  for  the  conditional  expectation  in  (11)  and  proved 
that  there  exists  a  one-to-one  transformation  of  ^  to  a  vector  (fn\  say,  of  independent  standard 
Gaussian  variables.  Since  the  transformation  is  a  bijection,  one  can  equivalently  use  in  (11)  C*n)  or 
C(n).  Consequently,  in  the  following,  we  shall  assume  that  is  transformed  to  a  vector  of  inde¬ 
pendent  standard  Gaussian  variables.  Generally,  for  random  times  Si,  the  conditional  expectation 
in  (11)  is  a  nonlinear  function  of  and  the  residual  process  An  is  dependent  on  however, 
given  values  of  An  is  Gaussian.  This  is  a  simple  consequence  of  the  definition  of  the  vector 
and  the  random  times  S\,...,Sn  (10),  since,  given  =  z^n\  Si(z^n~^)  =  Si,  i  =  l,...,n, 
are  fixed  time  points,  and  hence  the  conditional  distribution  of  A(s)  given  ^  =  z*n)  is  Gaussian 
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with  mean  zero  and  covariance 


(12)  r(s,  i ^n-V)  =  Cov{ A(s),  A(i)\A(t),  A '(t),  A{Sl ),...,  A(sn)) . 

ti 

Let  y  be  a  decomposable  process,  i.e.  y($)  =  F(s,X,  A (s)),  see  Definition  1,  where  A  is  a  vector 
■*  valued  Gaussian  process.  Assume  that  the  vector  of  random  times  5i,...,S„  in  (10)  is  given. 

and  let  ^  be  the  corresponding  vector  (10).  Denote  by  An($|z(n0,  zM  =  (r,ri,Zi,. . .  ,z„),  a 
zero-mean  Gaussian  process  with  covariance  function  (12).  Consider  the  process  y(.s|z(n0  defined 
as  follows 

(13)  y(s\zW)  =  F(s,X,  £7[A(^)|C<n)  =  z(n)]  +  An(s|z(">)). 

Since  AT,  A  are  independent  y(s|z(n0  =  ?/(s)l^n*  =  z^n\  (where  =  denotes  equility  in  distribution). 
In  order  to  use  (5)  with  y(s)  replaced  by  2/(s|z^),  let  yu(t)(s|z(n0  he  the  process  (4),  i.e. 

(14)  2/u«)(s|z<n))  =  F(s ,  '*!(<)),  S(A(S)|C<n)  =  z(n)]  +  An(S|z("))), 

where  p*)r(u(<.))  is  defined  by  (3).  e,  in  (14),  we  are  using  that  £[A(0|C^  =  =  r  and 

An(<|zW)  =  0. 

We  turn  now  to  the  second  formula  for  the  Tu-density. 

Theorem  4.  Let  y(s)  be  a  decomposable  process,  i.e.  y(s)  =  F(s,X,  A(s)),  as  in  Definition  1, 
where  A  is  a  k-dimensional  zero-mean  Gaussian  process.  Further,  assume  that  there  exist  k  -  1 
continuously  differentiable  mappings  from  R+  x  R  x  Rk  into  R,  (F\,...,Fk-\),  such  that  for  all 

x  G  R  and  all  s,  0  <  s  <  t,  except  possibly  for  a  finite  number  of  s, 

(15)  F(s,x,-)  =  ^F(5,z,  •),  Fi(s,x,  •),...  ,Fk-i(s,x,-)^J 

is  a  one-to-one  mapping  from  Rk  into  Rk.  Let  u(s),v(s),  u(s)  >  v(s),  be  continuously  differentiable 
barriers.  If  for  all  s,  0  <  s  <  t,  the  following  conditional  expectations  are  finite 

FOy'COs/'COlMO  =  «(0>j/(0  =  «(0]  <  +oo.  £(li/(02/'(5)i|y(0  =  «(0  >v(s)  =  v(s))  <  +00 - 
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then 


(16)  /r.w  =E[{m  -  A‘)fHy;S)\y(t)  =  «(!)) /„(■)(»«) 

-  J  E^ly'U)  -  u'(s))*(tj'(t)  -  •i,(l))  +  /(!/;S)/(o,J,(!/)|!/(l)  =  u(i), 

2/(5)  =  u(s)  f(u(t),u($))ds 

~  ^[(3/'(5)  “  v'(s))~ {y'(t)  -  u'(t))  +  r(y,S)I{0i3)(y)\y(t)  =  u(t), 

y{s)  =  u(s)  f(u(t),v(s))ds, 

where  S  =  Si,...,Sn  are  random  times  defined  by  (10)  and  f  is  a  joint  density  of  y(t),y(s).  The 
indicator  function  /(?,',  Si,...,  5n)  is  given  by  (€),  f(o,S)('j)  IS  the  indicator  function  defined  equal 
to  1  if  the  sample  path  ofy  does  not  cross  the  barriers  t;.  r  prior  to  time  s  and  equal  to  0  otherwise 
and  x+  =  maa:(0,;e),  x~  =  max(Q,-x). 

Proof:  Since  j/(5|z(n))  (13),  A'd  =  (r,r1,zi,...,z„),  is  decomposable,  then,  by  Theorem  3,  the 
conditional  density  of  Tu  given  =  z<ni  is  defined  by 

/TU<‘-'(,lzln>)  =  £f%.oW-M”>))(!>W”l)-«'('))+|s(i|2(n|)  =  '‘W]/»(^.i)W())- 

By  Lemma  2,  the  density  of  is  given  by 

(17)  /t„k<.>M*(”))  =  £[/(0, <> (Sf«< ej(-|*<")))I •  I**®’ )  -  «'W)+'/W<)|r), 


where  yu(t)  is  given  by  (14)  and  /(u(t)|r)  is  defined  by  (6).  We  are  also  using  that  sl(l)(¥nh,  = 
y'u(t)(t\^)  is  a  constant  variable  dependent  only  on  zf0)  =  (r,n). 

Now,  for  all  t,  i  -  1  Ar,(s;|z(rii)  =  0,  where  Si,...,sn  are  the  values  of  random  times 

Si,...  ,  Sn  given  ^  =  z(n),  (=  R(n+2lk,  and  hence  we  are  allowed  to  multiply  the  expectation 

in  (17)  by  the  indicator  -/(y«(i)(*!^n^ );  sj , . . . ,  5„)  (8),  which  is  a  function  of  z(n\  Further,  for  all 
z(n\  we  have 


£[/(o,t)(2/u)]  =  1  “  /  /ru(yj^)  :  /m,j(s)<!s, 

Jo 
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where  yu(-)  =  yu(t)(-\^n))  (14). 

Now,  since  An(-|z("))  is  a  ^-dimensional  Gaussian  process  and  (15)  holds,  one  can  prove,  see 
Theorem  2  in  [10],  that  the  absorptions  times  Tu(yu),  Tv(yu)  are  given  by  (1).  Finally,  by  combining 
(1),  (17-18),  multiplying  (17)  by  the  density  of  and  integrating  out  z^,  we  obtain  (16). 

3.3  Bounds  for  the  density  of  Tu.  Since  0  <  I(o,s){y)  <  1)  then,  for  n  >  1,  we  have  the  following 
upper  and  lower  bounds  for  the  density  of  Tu 

(19)  /+(t;n)  =E[(y'(t)  -  u'(t))+ I(y\S)\y(t)  =  u(t)]fy{t)(u(t)) 

fu  (t;n)  =fu{^n)-  J  E^[y'(s)  ~  u'(s))+  (y'(t)  -  iL'(t))+ I(ytS)\y(t)  =  u(t), 

y(s)  =  u(s)  f(u(t),u(s))ds 

-  J  ^[(j/'(3)-v/(s))-(2/'(0-u/(t))+J(j/;S)|j/(t)  =  ti(t), 

y(s)  =  v(s)  f(u(t),v(s))ds. 

Further,  for  n  -  0,  the  bounds  /+(t;  0),  f~(t\ 0)  are  obtained  by  replacing  in  (19)  the  indicator 
I(V,S)  by  1. 

We  turn  now  to  the  problem  of  choosing  the  vector  S.  Obviously,  for  any  fixed  t  and  n,  the  best 
choice  of  S  is  that  which  minimizes  the  upper  bound  /*(i;n),  or  the  difference  between  the  bounds 
/+(t;n)  -  /” (£;n),  see  (19).  However,  since  these  procedures  lead  to  complicated  optimization 
problems,  we  propose  a  simpler  recursive  procedure. 

We  begin  with  some  simple  properties  of  the  bounds  ft  if Z  (19)-  Assume  that  we  have  selected 
a  vector  of  random  times  S  =  (5i,...,5n),  and  let  be  a  random  vector  defined  by  (10).  In 
order  to  simplify  notation,  we  shall  denote  the  conditional  process  yu(t){' l2^)  in  (14)  by  yu.  Now, 
similarly  as  in  (17),  we  can  write  the  density  of  Tu  as  follows 

(20)  /r.(0  =  J  fTu^M^n))f^Mn))^ 

=  J  E[l{0il)(yu)]l(yu;si,...,sn)f(zw)f^)(z{n))d^n\ 
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where  yu(s)  =  J/u(f)(-s|^'1^),  s*  =  are  the  values  of  Si  for  ^'-1)  =  ^ and  f(z (0))  is 

defined  by 

/(^0)) =w(,  l(<i^0>)-«'(o)+/(>‘(i)i'-)- 

Further,  fan)  is  the  density  of^n).  Consequently,  the  bounding  problem  of  the  density  of  7’,,  is 
reduced  to  the  construction  of  an  upper  and  lower  bound  for  the  expectation 

(21)  E[I(o,t)(yu)]  =  P(v(s )  <  t/u(o('sMn))  <  u(5)  for  all  s.O  <  s  <  <)• 

Let  Pq,Pq  be  the  following  upper  and  lower-bounds  for  the  probability  (21) 

(22)  P0+(t\yu)=  1, 

Po(t-,Vu)  =1  -  /  £[(!/u(5)  ~  u'($))+\Vu{s)  =  u{s)]fyu{s)(u(s))ds 
Jo 

-  [  E[{y'u(s)  -  v\s)V\yu{s)  ==  v(s)]fyu{s)(v{s))ds, 

where  the  lower  bound  P0~  is  obtained  using  (1)  and  (18).  Now,  by  replacing  the  expectation  in  (20) 
by  the  upper  and  lower  bound  Po(t\yu),  Po(t;yu) ,  we  obtain  the  bounds  /+  (t\n),  /“(t;n)  (19). 
respectively,  generated  by  the  vector  =  (A(t),A'(t),A(Si),...,A(Sn))=  (^0),Cii- . - ,Cn) - 
In  order  to  obtain  more  accurate  bounds  f£(t\n+l),f~(t\n  +  l),  we  have  to  choose  an  additional 
random  time  5n+i,  0  <  Sn+ i  <  t,  which  is  a  function  of  tfnK  Note,  that  the  optimal  strategy  is 
to  select  the  whole  new  vector  ^n+1),  so  that  f+(t;n  +  1)  is  minimized.  However,  here  we  are 
restricting  ourselves  to  recursive  selection  procedures  of  5n+i,  i.e.  we  add  ^n+i  to  the  old  vector 
S. 

Now,  assume  that  we  have  selected  recursively  k  additional  random  times  Sn+i,-  •  . ,Sn+fc  and  let 
Cn+i,...)Cn+Jt  be  the  vector  ( A(5’n^.1 ),...,  A(5„+fc))  transformed  to  iid.  standard  Gaussian  vari¬ 
ables.  (The  selection  procedure  will  be  given  later  in  this  subsection.)  The  vector  (£n+1 , . . .  ,Cn+fc) 
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generates  new  bounds  for  the  probability  (21),  P*,Pjf,  say,  defined  as  follows 

(23)  P£(t;yu)=E[l(yu‘,Sn+x,...,Sn+k)\(<n)  =  *{n)j 

Pk(t',yu)=Pjk(t'<yu)~  I  E  (j/u(s)- «'(5))+/(j/u;5'„+i,...,5n+A:)|C(n)  =  *(n\ 

do 

yu(s)  =  n(s)  fyu(s){u(s))  ds 

-  f  E  [(^(s)  "  *'(*))' “/(ffu! 5n+i, .  •  • , S»+fc) |C<n)  =  *(n) , 

l lu(s)  =  v(s)  fyu(s)(v{s))ds. 

Once  again,  by  replacing  the  expectation  in  (20)  by  the  upper  and  lower  bound  P£ (t;  ?/u),  Pk 7 (<;  yu ) 
we  obtain  the  bounds  fZ{t;n  +  k ),  /“ (t\n  +  k)  (19),  respectively,  generated  by  the  vector  ^n+0. 
In  the  following  lemma,  we  give  a  recursive  formula  for  the  bounds  (23). 

Lemma  5.  The  upper  and  lower  bounds  (23)  P£  (j/u;t),  P(f  (yu\t),  k  >  0,  for  the  probability  (21), 
satisfy  the  following  recursive  formula 

(24)  P+(t;yu( o(-!z(n)))=  j  ^i(^2/.40(-lz(n+1))),/(^(0(-Mn+1));6"+i)/c»+.(z)(fe> 

Pk  (t; j/U(t)(-k(n)))= j  pk- 1  (*; 2/u(o(-l2(n+1))) -i{yu(t){¥n+1))\ 5n+i)/c„+, (*) 

where  sn+ i  =  5n+i(*^)  is  the  vaiue  of  5n+i  for^n)  =  and  z*n+1)  =  (z<n),z). 

Proof:  We  prove  the  lemma  only  for  the  upper  bound  P£ .  A  full  proof  is  only  notationally 
more  complicated.  By  additional  conditioning  on  Cn+i  in  (23)  and  using  (8), (13),  Pfc"(t;  2/u)  can  be 
written  as 

A+(‘i!M<>(-l*w))=/ . s„+*)l<<n)  =  *<n).C„+>  =*j/c.«to* 

=  J  E[l(,M<t(:\z<n*liy,S,lt, . S„+l..)|C<"+,)  =z<"+‘>] 

i[yu(t)(¥n+l));sn+i)fcn+l(z)  dz , 

showing  (24). 

Observe  that  in  the  recursive  definition  of  the  bounds  P*(<;  y.j),  Pk (t;  yu)  (24),  we  have  assumed 
that  the  random  times  5n+i, . . .  ,Sn+k  are  given.  Consequently,  in  order  to  use  (20)  and  (24)  to 
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calculate  /+(t;n  +  k),  f+(t\n  +  k),  we  have  to  define  a  recursive  procedure  to  choose  the  random 
time  Sn+i  as  a  function  of  In  addition,  since  yu(t)(>s|2vn^)  =  F(s,pt,r(u(t)),E[A(s)\^n'1  = 
z(n)]  +  An(s|z(n))),  where  the  only  random  component  is  a  zero-mean  vector  valued  Gaussian 
process  A„(s|^n^),  then,  once  the  random  time  S’„+i(^n^)  =  sn+i  is  chosen,  one  can  easily  obtain 
the  distribution  of  the  process  2/u(i)(sMn+1^)  by  calculating  An(s|zW)|An(sn+i|z(n))  =  z. 

We  turn  now  to  the  presentation  of  the  procedure  for  choosing  the  random  times  S  in  (19). 

P:.  Step  1:  choose  the  time  Si,  0  <  5i  <  t,  to  minimize  Pi’(t;yu(t)(,l^)),  and,  by  (23),  given  the 
values  z(°)  =  (r,ri),  choose  the  time  S\  to  minimize  P(v(s i)  <  2/u(t)(5iM°^)  <  u(si)). 

Step  n:  given  the  time  points  5i,...,5n-i  choose  the  time  Sn,  0  <  Sn  <  t,  to  minimize 
P+(i  Le ■  given  choose  the  time  sn  to  minimize  P(v(sn)  <  yu(t)(sn|^n-1*)  < 

“(«„))• 

Since,  for  any  z<n)  €  p(n+2)k,  n  >  0,  the  procedure  P  defines  j/u(t)(-|z(”))),P;f(f;T/u(f)(-|z<’l))). 
then,  using  (24),  we  can  recursively  evaluate  the  bounds  P*(t;z/u(f)(-|z(0))),  P,T(t;z/u(t)(-|z(0))). 
Hence,  the  upper  and  lower  bounds  /+(t;n),  /“(*;n)  are  defined  by 

(25)  /+((;«)  =  Jp+{t  i!/„(,)(|r|0|))/(2<0|)/(,.,(rl°>)iiz(0', 

/;(<;»)  =  J  PZ(> 

where  /(z^)  =  ( y'u{t\ z^)  -  u'(t))+ f(u(t) |r).  In  the  following  subsection  we  present  a  program 
BOUND,  which  evaluates  the  bounds  P+(t-,yu{l)(-\^)),  P~(t-,yu{t)(-\zi°'>))  in  (25),  for  the  special 
case  when  i/u(t)(-|z^)  is  a  zero-mean  Gaussian  process. 

3.4  Program  BOUND.  The  procedure  BOUND  evaluates  the  upper  and  lower  bound  for  the 
probability 

(26)  Prob  =  P(u(s)  >  A(s)  >  v(s)  for  all  s,0  <  s  <  t), 

where  A  is  a  continuously  differentiable  zero-mean  Gaussian  process.  We  begin  with  a  simple 
lemma. 
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Lemma  6.  Let  X  be  a  Gaussian  variable  with  mean  m  and  variance  <r2,  then 


JS[jr+]=«r-®(2), 


where  '$(x)  =  <f>(x)  +  x$(x),  <p  and  $  being  the  standardized  normal  density  and  distribution 
functions. 

In  all  calculations,  we  have  approximated  the  ^-distribution  by  Hermite  polynomials,  and  hence 
the  ^-function  is  very  accurately  approximated  by  an  explicit  function. 

We  turn  now  to  the  description  of  the  procedure  BOUND 

Input  variables: 

1:  n  -  number  of  iterations,  n  >  0, 

2:  t  -  fixed  time  point,  t  >  0, 

3:  u(s),v(s)  -  continuously  differentiable  barriers,  u(s)  >  u(s),  0  <  o  <  t, 

4:  u'(s),v'(s)  -  derivatives  of  the  barriers, 

5:  r(si,S2)  -  covariance  function  of  a  zero-mean  Gaussian  process  A,  0  <  si,S2  <  t, 

6:  ^(si,^) '  covariance  function  Cov(A'(si),  A(so))  =  gf^r(si,S2)>  0  <  si,S2  <  t, 

7:  o\ (s)  -  variance  of  the  derivative  A',  i.e.  <7f(s)  =  V ar(A'(s)),  0  <  s  <  t, 

Output  variables: 

Pn  '  upper  bound  for  the  probability  Prob  (26), 

P~  -  lower  bound  for  the  probability  Prob  (26). 

Algorithm : 

Since  the  bounds  P+,  P~  are  functions  of  u(-),v(-)  and  the  covariance  function  of  the  process 
A,  say,  we  express  it  in  notation  by  introducing  P+(u,  v,r&)  and  P~(u,v,Tb). 

If  n  =  0,  then  by  (22) 
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Po+  =  1, 


PQ  =  1-  /  E[(A'(s)  -  ti'(s))  +  [A(s)  =  u(s)]U{3){u(s))ds 
Jo 

~  [  ^[(A'(«)-'y/(s))~|A(5)  =  v(s)]fA(s)(v(s))ds. 
Jo 

Using  Lemma  6,  the  lower  bound  Pq  can  be  written  in  more  explicit  way,  i.e. 


0  Jo  o-(s)  \  *a(«)  0  ff(s) 


v'(s)-u(s)&(s)wv(sM 


a2(s) 


ds , 


where  the  ^-function  is  defined  in  Lemma  6  and  6(s)  =  ri(s,s)/r(s,s), 


£[A'(s)|A(s)  =  ?z(.s)]  =  -  u(s )  •  6(s), 

£[A'(.s)|A(,s)  =  u(s)]  =  -  v(s)  •  b(s), 

4M  =  r«r(A'(S)|AM)  =  <rj(»)  - 

Observe  that  the  integral  must  be  evaluated  numerically. 

If  n  >  1,  then,  by  the  procedure  P,  we  choose  $i,  0  <  si  <  t,  a  fixed  time  point  for  which  the 
probability  P(u(si )  >  A(si)  >  v($i))  is  minimized,  i.e. 

>  a(Si)  >  »(*,))  =  0™5,(*(^)  - 

and  by  recursion  (24)  we  have 

/■«(*  i) 

(27)  P*(«,v,rA)  =  /  Pn-iiu-xbiW^-xb^^^s^f^Wdx 

Jv(ai) 

N 

*'£,Pn-l{u-xMs)iv-xMs)’rA\M’i))fM'i)(Xi)hi’ 

i=l 

where  x,-,/i,'  are  suitable  nodes  and  weights,  respectively,  u(si)  >  x,-  >  v(si),  i  =  and 

(28)  x6i(s)  =  £{A(s)|A(si)  =  x]  =  x^~^y 
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The  same  recursion  can  be  given  for  the  lower  bound  P“,  by  replacing  in  (27)  ”+”  by 
Finally,  for  each  a we  evaluate  the  bounds  P*_x(«  -  Xi&i(s),t>  -  Xt&i(s),rA|A(Sl)),  P^i(n  - 
Xibi(s),v  -  Xibi(s),rA\&(Sl)),  using  the  procedure  BOUND  with  the  following  input  variables; 


(29) 


1:  n-  1, 

2:  t,  unchanged, 

3:  u(s)  -  Xibi(s),v(s)  -  x;6i(s),  where  &i(s)  is  given  by  (28),  0  <  s  <  t, 
4:  u'(s)  -  Xib[(s),v'(s)  -  Xib[(s),  where  b[(s)  is  given  by  (29),  0  <  s  <  t, 
5:  covariance  function  (30), 

0:  covariance  function  (31), 

7:  variance  (32), 

ri(si,s) 


b[(s)  = 


r(susi)’ 


(30) 


r&\A(3l)(s,t)  =  r(s,t)  - 


r(s,si)r(si,t) 
r($l}Sl)  ’ 


(31)  Co»(A'(.),A(l)|A(»,))=  r,(Sll)  - 


(32)  V'or(A'(«)|A(*,))=  <tj(js)  — 

The  output  variables  P+  and  P“  are  now  defined  by 

(33)  P*  —  Pjt-l  (U  ~  l(^)i v  ~  •gi^l(^)>rA|A(ji)l  /&(< 

i=l  '  ' 

Pn“  =Y,Pn-l[u  ~  XMS)^V  ~  xMs)^r^\Ms1))  fA(sx)(Xi)hi. 

Here,  we  assume  that  the  procedure  BOUND  is  programed  in  a  computer  language  which  allows 
recursive  functions,  e.g.  APL. 

Finally,  an  obvious  question  is  whether  one  should  use  in  (19)  fixed  times  si,...,s„  instead 
of  random  is  would  drasticly  reduces  the  number  of  times  one  have  to  evaluate 
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equations  (28-32).  However,  tests  disclosed  that  such  a  procedure  is  usually  much  slower.  The 
reason  is  that  the  probability  P(u(si)  >  A(si)  >  v(si))  >  P(u(S\)  >  A(Si)  >  i>(5i)),  and  hence 
we  have  to  increase  N,  number  of  nodes  a to  compute  (33). 

3.5  Concluding  remarks.  In  this  subsection  we  discuss  the  convergence  of  the  bounds  (19),  in 
the  case  when  S,  are  deterministic  and  dense  points  in  the  interval  [0,<].  Tiiis  does  not  prove  that 
the  bounds  obtained  using  the  procedure  P  converges  to  the  density  of  Tu.  However,  the  procedure 
evaluates  simultaneously  the  upper  and  lower  bounds  for  the  density  of  Tu ,  so  the  accuracy  and 
the  convergence  of  the  bounds  can  be  easily  checked.  In  addition,  in  examples  presented  in  the 
next  section,  the  bounds  based  on  the  random  times  S,  are  more  accurate  then  the  bounds  based 
on  deterministic  points  s,  and  the  algorithm  is  much  faster. 

Remark  7.  Under  assumptions  of  Theorem  4,  if  {s„}^  is  a  dense  subset  of  the  interval  [0,t], 
then  as  n  tends  to  infinity 

/(j/;sx,...,s„)  i/(0, <)(»/),  a .s., 

E  [(y'(0  ~  u'(0)+(y'(a)  -  u'(s))+/(y;  su...,  sn)\y(t)  =  u(t),  y(s )  =  «(«)]  1  0, 

E [(»'(0  ~  «'(0)+(2/'(a)  “  t/(s))~/(y; su..., Sn)|j/(0  =  «(t), y(s)  =  «(«)]  1  0- 

Consequently,  the  upper  and  lower  bounds  fZ(t;n),  /“(t;n),  (19),  with  Si  =  Si,  i  =  1,2,..., 
converge  monotonically  to  the  density  ofTu  as  n  goes  to  infinity.  (The  proof  is  similar  to  the  proof 
of  Theorem  11  in  [13].) 

Observe  that  Theorems  3  and  4  can  be  used  to  construct  many  different  types  of  upper  and  lower 
bounds.  For  example,  another  type  of  lower  bound  are  obtained  by  replacing  the  indicator  I(y;S) 
in  (16),  by  1,  and  overestimating  the  indicator  I(o,3)(y),  by  I(y;Si,...,Sn),  where  Si,...,Sn  can 

be  chosen  by  a  procedure  similar  to  P,  see  [13].  However,  an  important  property  of  (19),  which 

distinguishes  it  from  the  other  approaches,  is  that  the  same  points  S  are  used  in  both  upper  and 
lower  bounds  leading  to  a  more  efficient  algorithm. 
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4. Applications 


4.1  Slepian  model  process.  Let  y(t),  t  >  0,  be  a  stationary  zero-mean  ergodic  Gaussian  process 
with  covariance  function  r,  and  assume  that  its  sample  paths  are  a.s.  continuously  differentiable. 
A  sufficient  condition  [1]  for  this  is  that  the  process  is  separable  and  that 

r"(s)  =  A  -  o(|/o(7|s|ra), 

as  s  -*  0,  for  some  a  >  1.  Assume  that  the  process  y  is  normalized  so  that  Ao  =  Var(j/(0))  =  Ao  = 
Var(r/'(0))  =  1,  which  is  only  a  matter  of  scaling. 

In  following  subsections,  we  are  interested  in  the  ’’long  run”  properties  of  the  process  y  after 
downcrossings  of  the  level  u;  consequently  we  are  introducing  the  Slepian  model  process  £u  for  y 
after  a  downcrossing  of  the  level  u.  This  is  the  stochastic  process  £„(•)  which  is  distributed  as  the 
long  run  distribution  of  y(w,  tk  +  •),  when  tk  runs  over  all  w-downcrossings  of  t/(u>,  •) .  Mathematical 
details  about  Slepian  processes  and  long  run  probabilities  can  be  found  in  [5],  Ch.  10,  and  [6,  7]. 
We  now  give  a  simple  representation  of  the  Slepian  model  process  fu. 

Consider  a  zero-mean  Gaussian  process  A,  with  covariance  function 

(34)  Cov(A(s),  A(t))  =  Cov(y(s),  y(t)\y(0),y'(0))  =  r(s  -t)~  r(tf  -  r'(t )2, 

since  Var(y(0))  =  Var(?/'(0))  =  1,  and  let  R  be  a  standard  Rayleigh  variable,  with  mean 
independent  of  A.  Then  the  Slepian  model  process  £u  is  given  by 

(35)  £u  (s)  =  ur(s)  +  Rr'(s)  +  A(s). 

Obviously,  the  process  satisfies  the  assumptions  of  Theorem  4,  and  hence  we  can  use  the  bounds 
(19)  developed  in  Section  3.3. 

In  following  numerical  examples,  we  shall  use  the  process  y  with  covariance  functions  given  by 


4.2  The  distribution  of  wave-length  and  amplitude.  Assume  y(t)  is  a  zero-mean  Gaussian 
process,  which  described  water  elevation  at  a  fixed  point.  A  question  that  arises  in  oceanography 
is  that  of  the  ’’empirical”  or  ’’long  run"  distribution  of  the  zerocrossing  wave-length  and  amplitude 
T,  H.  By  this  we  mean  the  difference  between  the  time  of  the  zero  downcrossing  and  the  following 
zero  upcrossing  and  the  lowest  value  of  y  in  this  interval,  see  Figure  1.  Since  the  ’’long  run” 
properties  of  y  after  zero  downcrossings  are  described  by  Slepian  model  process  £o  (35),  we  have 

T,H 

For  notational  convenience,  we  drop  the  subscript  ”0”  in  and  in  the  following  the  Slepian  model 
process  after  zero  downcrossing  is  denoted  by  £. 


Figure  1.  Definition  of  zerocrossing  wavelength  and  amplitude  T,H. 


Obviously,  H($)  >  h  if  and  only  if  £(s),  s  >  0,  reaches  the  level  -h  before  it  crosses  the  level  0 
again.  (Observe  that  £(0)  =  0,  £'(0)  <  0.)  Thus  the  distribution  of  T,H  can  be  expressed  using 
the  densities  of  the  absorptions  times  of  f,  TU,TV,  with  u(s)  =  0,v(s)  =  -h,  respectively,  i.e. 


(37) 

(38) 


P(T  <t,H  <h)=  f  frAs)d3, 

Jo 

/•+ co  r+co 

P{H  <h)=  fTu(s)ds  =  1  -  /  fT.(s)ds. 
Jo  Jo 
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In  (38)  we  are  using  a  fact  that,  for  Gaussian  processes,  the  probability  that  f  stays  for  ever  between 
the  finite  barriers  is  zero,  i.e. 

r+oo 

/  (/t„(s)  +  /tu(s)]^  =  1. 

Jo 

Note  that  the  distribution  of  T  can  be  obtained  from  (37)  by  choosing  the  lower  barrier  v(s)  =  -oo. 

In  order  to  use  (38)  for  bounding  the  distribution  Fh,  we  have  to  approximate  the  infinite  region 
of  integration  by  some  bounded  interval.  For  many  processes  of  practical  interest,  there  exists  a 
positive  constant  To ,  such  that  fru(s)  ~  'n  (37),  for  all  s  >  To.  Now,  using  the  lower  bounds 
(19)  for  the  density  of  Tu,  u(s)  =  0,i>(s)  =  -oo,  respectively,  we  can  find  To  as  the  first  time  when 


fT° 

/  e, 

Jo 

for  some  small  e,  i.e.  P(T  >  T0)  <  e. 

We  turn  now  to  the  presentation  of  the  bounds  /+(s;  n),  f~(s ;  n )  for  the  density  of  Tu(f)  obtained 
using  the  procedure  P  of  Section  3.3.  The  bounds  for  the  T„(5)  density  can  be  derived  in  the  similar 
way. 

For  a  fixed  value  h ,  the  formula  (1)  for  the  density  of  Tu,  with  xi(s)  =  0,u(.s)  =  -h,  is  given  by 


(39)  /r.w  =  £[/(»,  .)(oe'w+|f;o  =  o]/{„)(o). 

The  formula  (39)  can  also  be  expressed  in  terms  of  y,  i.e. 

/t„(0  =  c  •  £[/(o,oW(°)V(0+|2K0)  =  0 ,2/(0  =  0]/y(0).y(o(°»0) 

/0  r+oo 

I  P( 0  >  y(s)  >  -h  for  all  5,0  <  s  <  i|y'(0)  =  z,y'(t)  =  zx, 

■oo  JO 

y( o)  =  0 ,y(t)  =  Q)f{z,zi)dzi  dz , 

where  c*"1  is  the  average  number  of  zero  downcrossings  per  unit  interval 
(40)  c  =  WorlrfO)  =  O)/,(o,(0)  =  (2 *)->  y/5, 
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by  the  celebrated  Rice  formula  [9],  and 


(41)  f{z,z  i)  =  \z\z\  f y1  (o)  ,y'  (t)  ,y(o)  i  ~i  1 0, 0),  z  <0  and  zx  >  0. 

In  the  following  we  assume  that  t,h,z,zi  are  fixed  values.  Many  of  formulas  will  depend  on 
t,h,z,zi,  however  for  notational  convenience  we  shall  not  always  write  this  dependence  explicitly. 
Let  m(s)  be  the  following  conditional  expectation 

m(s)  =  £[2/(s)|t/'(0)  =  z,y'{t)  =  zx,y(0)  =  0  ,y(t)  =  0], 

and  let  A  be  a  zero-mean  Gaussian  process  with  a  covariance  function  r(si,s2)  given  by 


r(si,s2)  =  Cov(y(si),  j/(so)|2/,(0),2/,(^),2/(0),  y(t)). 


It  means  that  m(-)  +  A(-)  =  y(-)|j/(0)  =  z,y'(t)  =  zuy(Q)  =  0 ,y(t)  =  0. 
As  before,  let  r  be  the  covariance  of  the  y  process.  Then  with 


/  1  -r"(t)  0  —r'(t) 

c(t)  _  -r"(0  1  r'(0  o 

M  j  “  0  r'(t)  1  r(t) 

\  - r\t )  0  r(t )  1 


c(s)  =  (-r'(s),  -r'(s  -  t),r(s),r(s  -  «)), 
the  mean  m(s )  and  a  covariance  function  f(si,S2)  are  given  by 

(*\ 

m(s)  =  c(s)C(i)-1  q  , 

\  0/ 

f(si,s2)  =  r($x  -  s2)  -  c(s1)C(t)-1c(s2):r- 
Using  the  process  A,  the  formula  for  the  density  of  Tu  can  be  written  as 

/0  y+co 

/  p(  -m(s)  >  A(s)  >  -h  -  m(s)  for  all  s, 0  <  s  <  t)  f(z,  zx )  dz\  dz , 

•co  JO 
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Figure  2.  Densities  of  Tu,  h  =  0.5, 1, 1.5, 2,  -foo,  for  the  covariance  (36). 
where  f{z,z{)  is  given  by  (41). 

Now,  upper  and  lower  bounds  (19)  /+(i,n),  f~(t,n),  respectively,  are  obtained,  using  the  pro¬ 
cedure  BOUND  of  Section  3.4,  by  over  and  under  estimating  the  probability  in  (42). 

We  turn  now  to  presentation  of  the  numerical  bounds  for  the  ’’long  run”  wave-length  and  ampli¬ 
tude  distribution  for  Gaussian  process  y  with  covariance  function  r  (36). 

Table  1  shows  bounds  /+ (t;  0), . . . ,  /+(<;  4),  /" (f ;  4), . . . ,  /" (t ;  1)  for  the  zerocrossing  wavelength 
density,  i.e.  u(s)  =  0,  v(s)  =  -oo,  for  the  process  y  with  covariance  (36).  We  can  see,  that  the 
upper  and  lower  bounds  are  almost  identical.  In  addition,  the  integral  of  the  lower  bound  /"(t;  4), 
over  the  interval  [0,12.5],  is  0.999,  indicating  that  only  0.1%  of  all  waves  are  longer  than  12.5. 
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0.00  2.00  4.00  6.00  3.00  10.00 


Figure  3.  Isolines  of  joint  density  of  wavelength  and  amplitude  T,H ,  for  the  covari¬ 
ance  (36). 


Table  1.  Bounds  /+(<;  0),  /+(«;  1),  /+(<; 2)  /+(<;3),  /+(*; 4),  /«  (<;4),  /-(*;3),  /„"(<;  2),  /“(<;  I), 


for  the  zerocrossing  wave-length  density,  i.e.  u(s )  =  0,t/(s)  =  -oo,  covariance  function  r  (36). 


t 

/u+(*;0) 

• 

• 

• 

/u+(*;4) 

/«(*; 4) 

• 

• 

/«(*;!) 

1 

0.123 

. 

. 

• 

. 

. 

0.123 

2 

0.352 

• 

• 

• 

• 

0.352 

3 

0.257 

• 

- 

•  _ 

. 

0.257 

4 

0.078 

0.073 

• 

• 

0.07 

0.072 

5 

0.114 

0.063 

0.062 

0.062 

0.0<*  L 

0.057 

6 

0.231 

0.071 

0.062 

0.056 

0.056 

0.056 

0.056 

7 

0.173 

0.039 

0.033 

0.027 

0.026 

0.026 

0.026 

0.025 

8 

0.109 

0.032 

0.021 

0.015 

0.014 

0.014 

0.014 

0.013 

9 

0.158 

0.059 

0.023 

0.017 

0.013 

0.013 

0.008 

0.000 

10 

0.205 

0.061 

0.019 

0.013 

0.009 

0.009 

0.005 

0.000 

11 

0.146 

0.040 

0.012 

0.009 

0.006 

0.004 

0.004 

0.001 

0.000 

12 

0.128 

0.043 

0.018 

0.008 

0.005 

0.002 

0.000 

0.000 

0.000 

In  order  to  bound  the  distribution  of  zerocrossing  amplitude  we  have  to  bound  density  of  Tu  for 
the  barriers  u(s)  =  0  and  v(s )  =  -h,  h  >  0.  For  a  fixed  t,  the  density  of  Tu  is  increasing  function  of 
h.  This  is  shown  on  Figure  2,  where  we  present  the  density  of  Tu  for  h  =  0.5, 1, 1.5, 2,  +co.  Finally, 


the  joint  density  of  wavelength  and  amplitude  T,  H  is  obtained  by  numerical  differentiation,  on  h. 


of  the  density  of  Tu,  see  Figure  3. 
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4.3  The  distribution  of  Rainflow  cycle  amplitude.  When  a  piece  of  metal  is  subjected  to  a 
periodically  varying  load  small  microscopic  inhomogeneities  can  develop  into  open  cracks,  leading 
to  tatigue  failure  after  a  random  amount  of  time.  The  distribution  of  fatigue  life  length  depends  on 
the  amplitudes  of  the  applied  ’’load  cycles”.  One  then  needs  a  rule  to  combine  the  damages  caused 
by  the  different  cycles.  The  most  commonly  used  damage  rule  is  due  to  Palmgren  &  Miner,  and 
postulates  that  the  total  damage  caused  by  a  stress  history  {St}  of  load  cycles  is 

*(0  , 

D(t)  =  ^n7' 

k=l  ^  5‘ 

where  the  sum  is  extended  over  all  cycles  completed  at  time  t  and  Na  is  the  median  cycle  life 
obtained  from  tests  with  constant  amplitude  s.  The  median  life  is  predicted  to  be  the  time  t  which 
makes  D(t )  greate.  an  or  equal  to  one.  In  most  situations,  the  median  cycle  life  N3  is  large, 
between  104  and  107,  and  therefore,  by  ergodicity  of  the  load  process,  the  fatigue  life  is  predicted 


A 


(43) 


c-  J  T^fs(s)(ls’ 


where  fs  is  the  density  of  the  crgodic  (long  run)  distribution  of  the  cycle  amplitude  Sk  and  c  is  a 
mean  number  of  a  cycles  counted  in  the  unit  interval  [0, 1). 

Dowling  [2]  has  studied  the  accuracy  of  the  predictors  of  the  fatigue  life  T  based  on  eight  most 
commonly  used  counting  methods,  and  finds  that  only  the  rainflow  cycle  (RFC)  counting  method 
leads  to  prediction  agreeing  with  actual  lifes. 

Due  to  the  great  importance  of  the  RFC-counting  method,  many  different  algorithms  have  been 
proposed  in  the  literature.  However,  most  of  them  have  a  complicated  ’’sequential”  structure  which 
makes  them  difficult  to  apply  when  their  statistical  properties  are  studied.  The  following  definition 
of  RFC-cycle,  given  in  [12],  is  more  convenient  for  statistical  analysis  of  long  run  properties  of  the 
RFC-cycles. 
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DEFINITION  8:.  Let  y(r),  -T  <  r  <  T,  be  a  load  function,  and  let  {<*},  with  <  t-i  <  0  < 
to  <  ti  <  ... ,  be  the  times  of  the  local  maxima  of  y(-).  For  a  local  maximum  at  time  ti,  let  tf  be 
the  time  for  the  first  upcrossing  after  t ,  of  the  level  y(t,)  (or  tf=T  if  no  such  upcrossing  exists 
for  t{  <  t  <  T),  and  let  t~  be  the  time  for  the  last  downcrossing  of  y(t{)  before  ti  (or  t~  =  -T  if 
no  such  downcrossing  exists  for  -T  <  r  <  t).  Let  the  lowest  minima  in  the  intervals  (t~ ,tx)  and 
( ti,tf )  occur  at  ti,  tr,  respectively,  and  let  t'  be  the  time  when  the  higher  of  the  minima  y(ti), 
y(tr)  occur,  i.e. 

{  tr  if  y{tl)  <  y{tr), 

*?  = 

(  ti  otherwise. 

The  RFC-count  attaches  to  a  maximum  at  time  t,  a  Rainflow  cycle  originating  at  ti,  defined  as 
a  pair  of  the  maximum  y(U)  and  minimum  y(t‘),  the  amplitude  of  the  cycle  is  given  by 


S,  —  y(t{)  y(ti), 


see  Figure  4.  Furthermore,  the  empirical  bivariate  distribution  of  the  local  maximum  M  =  y{t%) 
and  the  corresponding  RFC-minimum  m  =  y(tj)  is  defined  as 


y,T) 


i r{tj  6  [— r,  T1];  y(tj)  <  u,y(tj )  <  u} 

#{<;  e  i-T,T\} 


s 


Figure  4.  Definition  of  Rainflow  cycle. 
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For  some  functions  y(-),  the  empirical  distribution  F,\f>m(u,v\  y,T)  (44)  diverges  as  T  — >■  +co. 

However,  when  y  is  a  sample  path  of  an  ergodic  process  a  limit  of  FM>m(u,v]y,T),  as  T - hoo, 

exists  almost  surely  and  defines  a  bivariate  distribution  function  F^3m,  say.  Obviously,  since  the 
RFC-amplitude  S  =  M  -  m,  once,  knowing  the  ergodic  distribution  F^3m,  we  can  evaluate  the 
predictor  of  fatigue  life  (43).  In  the  following  we  shall  present  an  approximation  of  the  ergodic  RFC- 
distribution,  based  on  the  bounds  for  the  absorptions  times  TU,TV,  u(s)  =  u,v(s)  =  v,  presented 
in  previous  sections. 

Observe,  that  the  marginal  distribution  of  m  is  the  same  as  the  ergodic  distributions  of  the 
height  of  local  minima  and  for  Gaussian  processes  can  be  given  in  an  explicit  formula.  Hence,  the 
evaluation  of  the  ergodic  distribution  F^3m  is  equivalent  to  calculation  of  Per9(M  >  u,m  <  v),  i.e. 

(45)  Fffju,  v)  =  FZ'(v)  -  Per3(M  >  u,m  <  v). 


Now,  using  ergodic  properties  of  marked  point  processes,  see  Leadbetter  et  al.  [5],  Chapter  10,  the 
probability  Per9(M  >  u,  m  <  v)  is  given  by 

neroc  \.r  ^  ^  _.\  .  £(#{<i  €  [0,  1];  y(tf)  >  U,y(tm{)  <  t>}] 

P  (M>u,m<v) - WTilM}] - ’ 

where  t,  are  the  times  of  local  maxima,  see  Definition  8.  In  [14],  we  have  proved  that  the  mean 
!?[#{<;  G  [0, 1]; 3/(i*)  >  ui2/(^)  <  w}]  is  equal  to  the  mean  number  of  u-downcrossings,  by  y,  in  the 
interval  [0, 1],  which  are  followed  by  a  downcrossing  of  the  level  v  without  crossing  the  barrier  u  in 
between.  More  exactly,  for  a  fixed  u,v,  let  {s,},  s;  >  0,  be  a  sequence  of  downcrossings  of  the  level 
u,  then 


Per3(M  >  u  m<v)=c-  Sifl  e  ^0l  +  0  crosses  the  leveI  v  before  >  0)1 

=c-  P(Zu(t)  crosses  the  level  v  before  u,t  >  0), 
where  £u  is  the  Slepian  model  process  for  y  after  u-downcrossing  (35)  and  c  is  given  by 


mil 

£[#{( 


i 6 1°’  mi  _ 

.•  e  [o,i]>]  _  V  A4exp  ’ 
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where  A4  =  Var(y"( 0)),  see  (40). 

Finally,  by  the  definition  of  the  variable  Tu,  with  u(s)  =  u,v($)  =  v,  we  have 

r+co  f+<x> 

(46)  c~l -Per9(M  >  u,m  <  v)  =  /  frv(s)ds  =  1  -  /  fra{s)ds, 

Jo  Jo 

since,  for  Gaussian  processes,  the  probability  that  £u  stays  for  ever  between  the  finite  barriers  is 
zero.  As  in  the  previous  section,  we  have  to  approximate  the  infinite  region  of  integration  in  (46), 
by  some  finite  interval.  If  the  level  u  is  relatively  small,  e.g.  u  <  a,  <r2  =  Var(?/(0)),  or  the  levels 
u,v  are  close  to  each  other,  then  one  can  usually  find  a  constant  To,  such  that 


(47)  /  (/J"(-s;n)  +  /^(s;  n)]  ds  >  1  -  e, 

Jo 

for  some  small  positive  e,  where  /“,/“  are  lower  bounds  (19)  for  the  densities  of  TU)TV,  respectively. 
However,  in  the  case  of  the  high  positive  u  and  low  negative  v,  the  tails  of  densities  of  Tu  and  Tv 
become  very  long,  consequently  T0  is  also  large.  Hence,  in  order  to  find  suitable  T0,  we  have  to  use 
bounds  /“  (s;  n),/“(s;  n )  with  high  values  of  n,  which  causes  numerical  difficulties.  Consequently, 
in  the  following,  we  present  an  approximative  method  to  evaluate  the  probability  Perg(M  >  u,m< 
v )  in  the  case  when  levels  u,v  are  high  and  low,  respectively. 

As  in  the  previous  section,  we  assume  that  y  is  a  zero-mean  Gaussian  ergodic  process.  For  a 
positive  constant  To,  denote  by  Pu(To)  and  -P„(To)  the  following  truncated  integrals 


(48) 


Pu(To)=  [  °  fTu(s)ds, 
Jo 

Pv(To)=  [  °  fr.Wds. 
Jo 


By  the  definition  of  TV,TU  variables  the  integral  in  (46)  can  be  written  as  follows 


f  +  oo 

(49)  /  fr,  (s)  ds  =  PV(T0)  +  (1  -  PV(T0 )  -  PU{T0 ))  •  P(Zu{t)  crosses  the  level  v 

Jo 

before  u,<  >  To|fu(t)  stays  between  v,u  for  all  t,0  <  t  <  To). 

Obviously,  if  (47)  is  satisfied,  i.e.  (1  -  Pv(To)  -  PU(T b))  <  the  second  term  in  (49)  is  less  than  e 

and  can  be  disregarded. 
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It  is  well  known  that  for  Gaussian  processes,  see  Leadbetter  et.  al.  [5]  for  suitable  conditions,  the 
point  processes  of  downcrossings  of  levels  u  and  v  converges  to  independent  Poiss  ;n  processes,  as 

u  —>  +oo  and  v - co.  Furthermore,  by  (34),  if  the  covariance  function  r(t )  and  its  derivative  ?■'(<), 

of  the  process  y,  converge  to  zero  as  t  goes  to  infinity,  then  for  large  t,  we  have  £u(t)  x.  y(t ),  where 
«  denotes  approximative  equality  in  distribution.  Consequently,  we  propose  to  approximate  the 
conditional  probability  in  (49)  by  the  corresponding  probability  evaluated  for  independent  Poisson 
processes  with  the  same  crossing  intensities  as  fu,  i.e. 

(50)  Papp(T0)  =  f+°°  K(t)  •  exp"  A"(,)+A“ (J)  *  dty 

JTo 

where  the  intensities  A“,  Aj  are  given  by  Rice  formula,  i.e. 

K(t)  =m'umut)  =  v]/ut)(v) 
a i(t)  - E[t‘u(t)+\Ut )  =  «]/<.(«(»)• 

Now,  for  fixed  To,  by  replacing  the  conditional  probability  in  (49)  by  Papp(To),  we  obtain  an 
approximation  P(u,v,To),  say,  for  the  probability  Pera(M  >  u,m  <  v),  viz. 

(51)  P(a,  vtT0)  =  PV(T0)  +  (1  -  PV(T0)  -  PU(T0 ))  •  Papp(T0). 

Finally,  by  combining  (45)  and  (51)  we  obtain  an  approximation  of  the  joint  distribution  of  maxi¬ 
mum  M  and  the  RFC  minimum  m 

(52)  FM,m(W,To)  =  F%3(v)  -  P(u,v,T0). 

We  turn  now  to  the  numerical  example.  Let  the  covariance  function  of  the  process  y  be  given 
by  (36).  Figure  5  shows  the  approximation  P(u,v,To)  as  a  function  of  To,  for  u  =  2  and  v  =  0.5, 
0.,  -0.5,  -1.,  -1.5,  -2.,  -2.5,  -3..  We  can  see  that  P(u,v\To)  (51)  stabilizes  very  quickly,  indicating 
that  the  constant  To  can  be  chosen  as  low  as  5,  what  substantially  reduces  the  numerical  effort  to 
evaluate  the  probabilities  Pv(To),  Pu(To)  (48).  Figure  6  shows  the  level  curves  of  the  approximation 
(52),  To  =  10,  of  the  distribution  of  ( M,m )  covariance  (36). 
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Figure  5.  Approximations  P(u,v,To)  (51),  for  u  =  2  and  v  —  0.5,0. 3.,  as  a 
function  of  To,  for  covariance  function  (36). 


Figure  6.  Isolines  of  approximative  distribution  (52)  Ft\ftm(v.,v]To),  To  =  10,  of 
maximum  M  and  the  RFC  minimum  m,  for  covariance  function  (36). 

Recently,  Ford  [4]  and  Nielsen  [7]  have  proposed  approximations  of  the  joint  distribution  of 


Figure  7.  Isolines  of  approximative  distribution  (52)  FM,m{u,v;T0),  To  =  0,  of 
maximum  M  and  the  RFC  minimum  m,  for  covariance  function  (36). 

( M,m )  equivalent  to  F^im(u,v}Q).  Since  FW,m(u,  v;  0)  approximation  is  based  on  the  assumption 
that  u-upcrossings  and  v-downcrossings  of  Slepian  process  are  independent  Poisson  point  pro¬ 
cesses,  this  approximation  can  be  accurate  only  for  high  positive  u  and  low  negative  v,  see  Figure 
7. 
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